clear
est clear


cap mat define migrationfigure=J(105,10,.)
		local row=1

	use "$input/full_geoid_pair_panel.dta"
	
	
			
	replace migration=migration/(tot_pop/100000)

	//Make sure only include states that have a storm in at least one county at least once in our sample years
	bysort state_fips: egen max=max(storm)
		drop if max!=1
		
	//Create future pooled migration - total migrants in the subsequent five years
	sort geoid_pair move_year
	g post_storm=1 if storm==1|storm[_n-1]==1&geoid_pair==geoid_pair[_n-1]|storm[_n-2]==1&geoid_pair==geoid_pair[_n-2]|storm[_n-3]==1&geoid_pair==geoid_pair[_n-3]|storm[_n-4]==1&geoid_pair==geoid_pair[_n-4]|storm[_n-5]==1&geoid_pair==geoid_pair[_n-5]
	mvencode post_storm, mv(0)
	
	//Collapse to total migration from a county (not looking at paired county to county migration, just all migration out of a county)
	collapse (sum) migration (mean) share_approved totalapprovedihpamount,by(from_geoid move_year storm post_storm state_fips)
	
	//Create IHS transformations
	ihstrans(migration)
	
	
	//Create lags of storms (was there a storm in t-1 years)

			
	label var post_storm "1(0-5 yrs post storm)"
	label var storm "1(Storm year)"
	sort from_geoid move_year
	
	replace totalap=totalap/1000000
	g highdamage=(totalapp>10 & totalapp<20) if totalapp!=.
		replace highdamage=. if totalapp>20
	mvencode highdamage, mv(0) o
		sort from_geoid move_year
		foreach i in 1 2 3 4 5{
				g highdamage_L`i'=highdamage[_n-`i'] if from_geoid==from_geoid[_n-`i'] 
		}	
	foreach var in L1 L2 L3 L4 L5{
		label var highdamage_`var' "1(Storm year), `var'"
	}
		
		reghdfe ihs_migration i.highdamage, absorb(i.from_geoid  i.move_year i.state_fips#c.move_year) vce(cluster state_fips)
		
			eststo t1c1
			cap mat def migrationfigure[`row',1]=_b[1.highdamage]
			cap mat def migrationfigure[`row',2]=_se[1.highdamage]
			cap mat def migrationfigure[`row',3]=15
			local row=`row'+1
		
		drop highdamage*
		g highdamage=(totalapp>20 & totalapp<40) if totalapp!=.
			replace highdamage=. if totalapp>40	
		mvencode highdamage, mv(0) o
		sort from_geoid move_year
		foreach i in 1 2 3 4 5{
				g highdamage_L`i'=highdamage[_n-`i'] if from_geoid==from_geoid[_n-`i'] 
		}	
		foreach var in L1 L2 L3 L4 L5{
			label var highdamage_`var' "1(Storm year), `i'"
		}
		reghdfe ihs_migration i.highdamage, absorb(i.from_geoid  i.move_year i.state_fips#c.move_year) vce(cluster state_fips)
	
			eststo t1c2
			cap mat def migrationfigure[`row',1]=_b[1.highdamage]
			cap mat def migrationfigure[`row',2]=_se[1.highdamage]
			cap mat def migrationfigure[`row',3]=30
			local row=`row'+1
		
		drop highdamage*
		g highdamage=(totalapp>40 & totalapp<80) if totalapp!=.
			replace highdamage=. if totalapp>800
		mvencode highdamage, mv(0) o
		sort from_geoid move_year
		foreach i in 1 2 3 4 5{
				g highdamage_L`i'=highdamage[_n-`i'] if from_geoid==from_geoid[_n-`i'] 
		}	
		foreach var in L1 L2 L3 L4 L5{
			label var highdamage_`var' "1(Storm year), `i'"
		}
		reghdfe ihs_migration i.highdamage, absorb(i.from_geoid  i.move_year i.state_fips#c.move_year) vce(cluster state_fips)
			

		eststo t1c3
		cap mat def migrationfigure[`row',1]=_b[1.highdamage]
			cap mat def migrationfigure[`row',2]=_se[1.highdamage]
			cap mat def migrationfigure[`row',3]=60
			local row=`row'+1		
		
		drop highdamage*
		g highdamage=(totalapp>80 & totalapp<160) if totalapp!=.
			replace highdamage=. if totalapp>160
		mvencode highdamage, mv(0) o
		sort from_geoid move_year
		foreach i in 1 2 3 4 5{
				g highdamage_L`i'=highdamage[_n-`i'] if from_geoid==from_geoid[_n-`i'] 
		}	
		foreach var in L1 L2 L3 L4 L5{
			label var highdamage_`var' "1(Storm year), `i'"
		}
		reghdfe ihs_migration i.highdamage, absorb(i.from_geoid  i.move_year i.state_fips#c.move_year) vce(cluster state_fips)

		eststo t1c4
		cap mat def migrationfigure[`row',1]=_b[1.highdamage]
			cap mat def migrationfigure[`row',2]=_se[1.highdamage]
			cap mat def migrationfigure[`row',3]=120
			local row=`row'+1	
		
		
		drop highdamage*
		g highdamage=(totalapp>160) if totalapp!=.
		mvencode highdamage, mv(0) o
		sort from_geoid move_year
		foreach i in 1 2 3 4 5{
				g highdamage_L`i'=highdamage[_n-`i'] if from_geoid==from_geoid[_n-`i'] 
		}	
		foreach var in L1 L2 L3 L4 L5{
			label var highdamage_`var' "1(Storm year), `i'"
		}
		reghdfe ihs_migration i.highdamage, absorb(i.from_geoid  i.move_year i.state_fips#c.move_year) vce(cluster state_fips)

		eststo t1c5
		cap mat def migrationfigure[`row',1]=_b[1.highdamage]
			cap mat def migrationfigure[`row',2]=_se[1.highdamage]
			cap mat def migrationfigure[`row',3]=200
			local row=`row'+1
			

label var highdamage "1(Storm year)"				
				
	preserve
		clear
			
		svmat2 migrationfigure

		drop if migrationfigure1==.
		
		
		
		g HI95_mean=(1.96*migrationfigure2)+migrationfigure1
		g LI95_mean=migrationfigure1-(1.96*migrationfigure2)

		
		
		set scheme plotplain
		
		rename migrationfigure1 mean_beta
		
		
		
		replace mean_beta=mean_beta*100 
		replace HI95_mean=HI95_mean*100 
		replace LI95_mean=LI95_mean*100
		
		
		reg mean_beta migrationfigure3
		
		twoway 		(lfit mean_beta migrationfigure3  , color(navy) lp(dash)) ///
					(scatter mean_beta migrationfigure3  , color("252 89 16%50")) ///
					(rspike HI95_mean LI95_mean migrationfigure3 , lc(gs8)  ///
					yti("Percentage change in out migration",size(vlarge)) ylabel(, labsize(large)) xlabel(, labsize(large))  ysize(1.5)  xsize(5) xtitle("Compensated damage (MM)", size(vlarge)) legend(off) )
					
					graph export "$figures/fig_si4.png", as(png) replace
					
					
	restore		
	
	
	
	
	
			